rm(list=ls())

library(terra)
library(tmap)
library(tmaptools)
library(sf)
library(leaflet)
library(haven)
library(rgdal)
library(maptools)
library(sp)
library(raster)

wbluescale <- c("#ffffff", "#08519c")
worangescale <- c("#ffffff", "#fc9272")
wgreenscale <- c("#ffffff", "#006d2c")
wgrayscale <- c("#ffffff", "#AAAAAA")
wpinkscale <- c("#ffffff", "#FF66B2")
wpurpscale <- c("#ffffff", "#9e9ac8")
wredscale <- c("#ffffff", "#de2d26")
bwhitescale6 <- c("#08519c", "#3182bd", "#6baed6", "#9ecae1", "#c6dbef", "#ffffff")
wbluescale6 <- c("#ffffff", "#c6dbef", "#9ecae1", "#6baed6", "#3182bd", "#08519c")

setwd("/Users/cb2257/Desktop/APSR Replication/")

ian <- st_read("Shapefiles/hurricane_sample/hurricane_sample.shp")
st_is_valid(ian, reason = TRUE)
ian <- st_make_valid(ian)

south <- st_read("Shapefiles/southern_states/southern_states.shp")
st_is_valid(south, reason = TRUE)
south <- st_make_valid(south)

eyepath <- st_read("Shapefiles/ian_eyepath/ian_eyepath.shp")
st_is_valid(eyepath, reason = TRUE)
eyepath <- st_make_valid(eyepath)

ida <- st_read("Shapefiles/ida_eyepath/ida_eyepath.shp")
st_is_valid(ida, reason = TRUE)
ida <- st_make_valid(ida)

nicole <- st_read("Shapefiles/nicole_eyepath/nicole_eyepath.shp")
st_is_valid(nicole, reason = TRUE)
nicole <- st_make_valid(nicole)

county <- read_dta(file="Data/Survey/map_descriptives.dta")

joined <- merge(ian, county, by.x="GEOID", by.y="geoid")

# Map of Respondents -- Figure 1 in Main Paper

respondents <- tm_shape(joined)+
  tm_borders()+
  tm_fill("response", title = "# of Respondents", palette=wgreenscale, style="fixed",
  breaks=c(0, 1, 4, 11, 27, 280), interval.closure="left", 
  labels=c("0", "1-3", "4-10", "11-26", "27+"))+
  tm_layout(" ", title.size=1, frame = F)

respeye<-respondents+tm_shape(eyepath)+tm_lines(col = "red", lwd = 2, lty = "dashed")

combined<-tm_shape(south)+tm_fill(col = "gray83")+tm_borders(col = "black", lwd = 1, lty = "solid")+respeye
tmap_save(combined, "Figures/respondent_map.png")

# Map of Respondents2 -- Figure A-12 in Main Appendix

respondents2 <- tm_shape(joined)+
  tm_borders()+
  tm_fill("response2", title = "# of Respondents", palette=wgreenscale, style="fixed",
          breaks=c(0, 1, 4, 11, 27, 58), interval.closure="left", 
          labels=c("0", "1-3", "4-10", "11-26", "27+"))+
  tm_layout(" ", title.size=1, frame = F)

respeye<-respondents2+tm_shape(eyepath)+tm_lines(col = "red", lwd = 2, lty = "dashed")

combined<-tm_shape(south)+tm_fill(col = "gray83")+tm_borders(col = "black", lwd = 1, lty = "solid")+respeye
tmap_save(combined, "Figures/respondent2_map.png")

# Map of Windswath -- Figure 2c in Main Paper

windswath <- tm_shape(joined)+
  tm_borders()+
  tm_fill("ian_windswath", title = "Windswath", palette=wbluescale, style="fixed",
  breaks=c(0, 1, 2, 3, 4), interval.closure="left", 
  labels=c("Sub-Cyclonic Winds (0-33 Knots)", "Tropical Storm Winds (34+ Knots)", "Violent Storm Winds (50+ Knots)", "Hurricane Winds (64+ Knots)"))+
  tm_layout(" ", title.size=1, frame = F)

windeye<-windswath+tm_shape(eyepath)+tm_lines(col = "red", lwd = 2, lty = "dashed")

combined<-tm_shape(south)+tm_fill(col = "gray83")+tm_borders(col = "black", lwd = 1, lty = "solid")+windeye
tmap_save(combined, "Figures/windswath_map.png")

# Map of FEMA Aid -- Figure 2f in Main Paper

fema <- tm_shape(joined)+
  tm_borders()+
  tm_fill("fema_aid", title = "FEMA Assistance", palette=wbluescale, style="fixed",
  breaks=c(0, 1, 2, 3, 4, 5), interval.closure="left", 
  labels=c("None", "Public (Category B)", "Public (Category A-G)", "Individual and Public (Category B)", "Individual and Public (Category A-G)"))+
  tm_layout(" ", title.size=1, frame = F)

femaeye<-fema+tm_shape(eyepath)+tm_lines(col = "red", lwd = 2, lty = "dashed")

combined<-tm_shape(south)+tm_fill(col = "gray83")+tm_borders(col = "black", lwd = 1, lty = "solid")+femaeye
tmap_save(combined, "Figures/fema_map.png")

# Map of Distance -- Figure 2d in Main Paper

distance <- tm_shape(joined)+
  tm_borders()+
  tm_fill("ian_distance_mi", title = "Distance to Eyepath", palette=bwhitescale6, style="fixed",
  breaks=c(0, 77.05172, 190.0623, 895.7505, 1154.872, 1393.435, 2000), interval.closure="left", 
  labels=c("< 10th Percentile", "10-25th Percentile", "25-50th Percentile", "50-75th Percentile", "75-90th Percentile", "> 90th Percentile"))+
  tm_layout(" ", title.size=1, frame = F)

distanceye<-distance+tm_shape(eyepath)+tm_lines(col = "red", lwd = 2, lty = "dashed")

combined<-tm_shape(south)+tm_fill(col = "gray83")+tm_borders(col = "black", lwd = 1, lty = "solid")+distanceye
tmap_save(combined, "Figures/distance_map.png")

# Map of Evacuations -- Not presented in paper

evacuation <- tm_shape(joined)+
  tm_borders()+
  tm_fill("ian_evacuation", title = "Evacuation Order", palette=wbluescale, style="fixed",
          breaks=c(0, 1, 2), interval.closure="left", 
          labels=c("No", "Yes"))+
  tm_layout(" ", title.size=1, frame = F)

evaceye<-evacuation+tm_shape(eyepath)+tm_lines(col = "red", lwd = 2, lty = "dashed")

combined<-tm_shape(south)+tm_fill(col = "gray83")+tm_borders(col = "black", lwd = 1, lty = "solid")+evaceye
tmap_save(combined, "Figures/evacuation_map.png")

# Map of Deaths -- Not presented in paper

deaths <- tm_shape(joined)+
  tm_borders()+
  tm_fill("i_ian_deaths", title = "Fatalities", palette=wbluescale, style="fixed",
          breaks=c(0, 1, 2), interval.closure="left", 
          labels=c("No", "Yes"))+
  tm_layout(" ", title.size=1, frame = F)

deathseye<-deaths+tm_shape(eyepath)+tm_lines(col = "red", lwd = 2, lty = "dashed")

combined<-tm_shape(south)+tm_fill(col = "gray83")+tm_borders(col = "black", lwd = 1, lty = "solid")+deathseye
tmap_save(combined, "Figures/i_deaths_map.png")

# Map of Deaths -- Figure 2e in Main Paper

deaths <- tm_shape(joined)+
  tm_borders()+
  tm_fill("ian_deaths", title = "# of Fatalities", palette=wbluescale, style="fixed",
          breaks=c(0, 1, 2, 6, 10, 73), interval.closure="left", 
          labels=c("0", "1", "2-5", "6-9", "10+"))+
  tm_layout(" ", title.size=1, frame = F)

deathseye<-deaths+tm_shape(eyepath)+tm_lines(col = "red", lwd = 2, lty = "dashed")

combined<-tm_shape(south)+tm_fill(col = "gray83")+tm_borders(col = "black", lwd = 1, lty = "solid")+deathseye
tmap_save(combined, "Figures/deaths_map.png")

# Map of Storm Surge -- Figure 2b in Main Paper

surge <- tm_shape(joined)+
  tm_borders()+
  tm_fill("ian_surge", title = "Storm Surge", palette=wbluescale, style="fixed",
  breaks=c(0, 1, 2, 3, 4, 5), interval.closure="left", 
  labels=c("None", "1-3 Feet", "3-6 Feet", "6-9 Feet", "9+ Feet"))+
  tm_layout(" ", title.size=1, frame = F)

surgeeye<-surge+tm_shape(eyepath)+tm_lines(col = "red", lwd = 2, lty = "dashed")

combined<-tm_shape(south)+tm_fill(col = "gray83")+tm_borders(col = "black", lwd = 1, lty = "solid")+surgeeye
tmap_save(combined, "Figures/surge_map.png")

# Map of Index -- Figure 2a in Main Paper

index <- tm_shape(joined)+
  tm_borders()+
  tm_fill("hurricane_pctile", title = "Hurricane Exposure Index", palette=wbluescale, style="fixed",
  breaks=c(0, 1, 2, 3, 4, 5, 6), interval.closure="left", 
  labels=c("< 10th Percentile", "10-25th Percentile", "25-50th Percentile", "50-75th Percentile", "75-90th Percentile", "> 90th Percentile"))+
  tm_layout(" ", title.size=1, frame = F)

indexeye<-index+tm_shape(eyepath)+tm_lines(col = "red", lwd = 2, lty = "dashed")

combined<-tm_shape(south)+tm_fill(col = "gray83")+tm_borders(col = "black", lwd = 1, lty = "solid")+indexeye
tmap_save(combined, "Figures/index_map.png")

# Map of Binary Index -- Figure A-5 in Main Appendix

binary <- tm_shape(joined)+
  tm_borders()+
  tm_fill("i_index", title = "Hurricane Exposure", palette=wbluescale, style="fixed",
          breaks=c(0, 1, 2), interval.closure="left", 
          labels=c("No", "Yes"))+
  tm_layout(" ", title.size=1, frame = F)

binaryeye<-binary+tm_shape(eyepath)+tm_lines(col = "red", lwd = 2, lty = "dashed")

combined<-tm_shape(south)+tm_fill(col = "gray83")+tm_borders(col = "black", lwd = 1, lty = "solid")+binaryeye
tmap_save(combined, "Figures/i_index_map.png")

# Map of Ida -- Figure A-9 in Main Appendix

idaindex <- tm_shape(joined)+
  tm_borders()+
  tm_fill("ida_pctile", title = "Hurricane Exposure Index", palette=wbluescale, style="fixed",
          breaks=c(0, 1, 2, 3, 4, 5, 6), interval.closure="left", 
          labels=c("< 10th Percentile", "10-25th Percentile", "25-50th Percentile", "50-75th Percentile", "75-90th Percentile", "> 90th Percentile"))+
  tm_layout(" ", title.size=1, frame = F)

idaeye<-idaindex+tm_shape(ida)+tm_lines(col = "red", lwd = 2, lty = "dashed")

combined<-tm_shape(south)+tm_fill(col = "gray83")+tm_borders(col = "black", lwd = 1, lty = "solid")+idaeye
tmap_save(combined, "Figures/ida_map.png")

# Map of Nicole  -- Figure SI-3 in Additional Supplement

nicoleindex <- tm_shape(joined)+
  tm_borders()+
  tm_fill("nicole_pctile", title = "Hurricane Exposure Index", palette=wbluescale, style="fixed",
          breaks=c(0, 1, 2, 3, 4, 5, 6), interval.closure="left", 
          labels=c("< 10th Percentile", "10-25th Percentile", "25-50th Percentile", "50-75th Percentile", "75-90th Percentile", "> 90th Percentile"))+
  tm_layout(" ", title.size=1, frame = F)

nicoleeye<-nicoleindex+tm_shape(nicole)+tm_lines(col = "red", lwd = 2, lty = "dashed")

combined<-tm_shape(south)+tm_fill(col = "gray83")+tm_borders(col = "black", lwd = 1, lty = "solid")+nicoleeye
tmap_save(combined, "Figures/nicole_map.png")

# Map of 2020 Republican Voteshare -- Not presented in paper

repvote <- tm_shape(joined)+
  tm_borders()+
  tm_fill("repvote20", title = "2020 Republican Voteshare", palette=wredscale, style="fixed",
  breaks=c(0, .2, .4, .6, .8, 1), interval.closure="left", 
  labels=c("< 20%", "20-39%", "40-59%", "60-79%", "> 80%"))+
  tm_layout(" ", title.size=1, frame = F)

combined<-tm_shape(south)+tm_fill(col = "gray83")+tm_borders(col = "black", lwd = 1, lty = "solid")+repvote
tmap_save(combined, "Figures/repvote_map.png")

repwon <- tm_shape(joined)+
  tm_borders()+
  tm_fill("repvote20", title = "2020 Presidential Election", palette=wredscale, style="fixed",
  breaks=c(0, .5, 1), interval.closure="left", 
  labels=c("Biden Won", "Trump Won"))+
  tm_layout(" ", title.size=1, frame = F)

combined<-tm_shape(south)+tm_fill(col = "gray83")+tm_borders(col = "black", lwd = 1, lty = "solid")+repwon
tmap_save(combined, "Figures/repwon_map.png")

# Map of Relief Aid -- Not presented in paper

aid <- tm_shape(joined)+
  tm_borders()+
  tm_fill("ihpAmount", title = "FEMA Individual \n Assistance per Applicant", palette=wpinkscale, style="fixed",
          breaks=c(0, .00000001, 495, 766, 1700, 2453, 2958), interval.closure="left", 
          labels=c("< 10th Percentile", "10-25th Percentile", "25-50th Percentile", "50-75th Percentile", "75-90th Percentile", "> 90th Percentile"))+
  tm_layout(" ", title.size=1, frame = F)

combined<-tm_shape(south)+tm_fill(col = "gray83")+tm_borders(col = "black", lwd = 1, lty = "solid")+aid
tmap_save(combined, "Figures/femaind_map.png")

aid <- tm_shape(joined)+
  tm_borders()+
  tm_fill("numberobligatedprojects", title = "FEMA Public \n Assistance Projects", palette=wpinkscale, style="fixed",
          breaks=c(0, 1, 2, 8, 12, 21, 84), interval.closure="left", 
          labels=c("< 10th Percentile", "10-25th Percentile", "25-50th Percentile", "50-75th Percentile", "75-90th Percentile", "> 90th Percentile"))+
  tm_layout(" ", title.size=1, frame = F)

combined<-tm_shape(south)+tm_fill(col = "gray83")+tm_borders(col = "black", lwd = 1, lty = "solid")+aid
tmap_save(combined, "Figures/femaproj_map.png")

aid <- tm_shape(joined)+
  tm_borders()+
  tm_fill("obligated", title = "FEMA Public \n Assistance per Project", palette=wpinkscale, style="fixed",
          breaks=c(0, .00000001, 15077, 47303, 150685, 639896, 3729566), interval.closure="left", 
          labels=c("< 10th Percentile", "10-25th Percentile", "25-50th Percentile", "50-75th Percentile", "75-90th Percentile", "> 90th Percentile"))+
  tm_layout(" ", title.size=1, frame = F)

combined<-tm_shape(south)+tm_fill(col = "gray83")+tm_borders(col = "black", lwd = 1, lty = "solid")+aid
tmap_save(combined, "Figures/femapub_map.png")
